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Abstract 

We perform Monte Carlo simulations of a super symmetric matrix model, which is 
obtained by dimensional reduction of 4D SU(iV) super Yang-Mills theory. The model 
can be considered as a four-dimensional counterpart of the IIB matrix model. We 
extract the space-time structure represented by the eigenvalues of bosonic matrices. 
In particular we compare the large N behavior of the space-time extent with the 
result obtained from a low energy effective theory. We measure various Wilson loop 
correlators which represent string amplitudes and we observe a nontrivial universal 
scaling in N. We also observe that the Eguchi-Kawai equivalence to ordinary gauge 
theory does hold at least within a finite range of scale. Comparison with the results 
for the bosonic case clarifies the role of supersymmetry in the large N dynamics. It 
does affect the multi-point correlators qualitatively, but the Eguchi-Kawai equivalence 
is observed even in the bosonic case. 
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1 Introduction 



A recent excitement in string theory is that we finally arrive at concrete proposals for non- 
perturbative definitions of superstring theory [l], @) H 3- m particular, Matrix Theory [|TJ 
and the IIB matrix model ||, which are candidates for a nonperturbative definition of M- 
theory and type IIB superstring theory respectively, have attracted considerable interest 
(for reviews, see Refs. j^, W). The proposed formulations take the form of large N reduced 
models J7j, which can be obtained by dimensional reduction of large N gauge theories; a 
reduction to one dimension for M-theory, and to zero dimension (one point) for type IIB 
superstring theory. These proposals are supported by some evidences such as the similar- 
ity of the Hamiltonian (or the action) to that of membranes or strings, the appearance of 
soliton-type objects known as D-branes with consistent interactions, and the consistency 
with string dualities upon compactification. For the IIB matrix model, even an attempt to 
establish a direct connection to perturbative string theory has been made by deriving the 
light-cone string field Hamiltonian from loop equations of the model ||. This attempt was 
indeed successful, albeit with the aid of symmetry and power-counting arguments. Further 
variants of that model have been proposed in Refs. || . 

As another approach to analyze these proposals we can investigate the dynamical prop- 
erties of large N reduced models of this kind, and verify if they really have the potential to 
describe nonperturbative string theory. In Ref. [TD| a two-dimensional reduced model with 
unitary matrices has been studied in this context. There, a large N limit, which differs from 
the planar limit (or 't Hooft limit), has been discovered numerically [J A Hermitian matrix 
model obtained by simply omitting the fermions in the IIB matrix model, and its general- 
izations to arbitrary dimensions larger than two, have been studied in Refs. [ p~2| , [13] . In Ref. 
T3| , Monte Carlo simulations up to N = 256 have been reported and analytical methods 



such as perturbation theory, Schwinger- Dyson equations and 1/D expansions have been ap- 
plied, providing a comprehensive understanding of the large N dynamics of that model. A 
new type of Monte Carlo technique was used to extract the value of the partition function 



14] , |12 |. This technique has been further applied to extract the asymptotic behavior of the 



eigenvalue distribution for large eigenvalues |T5|, |I6| . 

In the present paper, we make a first attempt to extract the large N dynamics of a super- 
symmetric large N reduced model obtained by dimensional reduction of 4D supersymmetric 



1 This non-planar large N limit has recently been re-interpreted as a continuum limit of non-commutative 
gauge theory Q. 
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Yang-Mills theory. The model can be regarded as a 4D counterpart of the IIB matrix model. 
The bosonic model is well understood [|l^, [13|], but the inclusion of fermions makes the system 
far more complicated. An attempt to study the model analytically maps it onto a soluble 
system |T7|], but the relevance to the original model is unclear due to a nontrivial change 
of variables in an analytic continuation. Here we take a direct approach, based on Monte 
Carlo simulations. Fermions are completely included by the use of the so-called Hybrid-R 



algorithm [I8|, which is one of the standard methods in QCD simulations with dynamical 
quarks. 

One of the features that makes the IIB matrix model most attractive as a nonperturbative 
definition of string theory is that space-time is dynamically generated as the eigenvalue 
distribution of the bosonic matrices 0, I5[ p0[| . In Ref. |1| a low energy effective theory 



of the model is constructed, where the authors discuss some possible mechanisms that may 
induce a collapse of the eigenvalue distribution to a four-dimensional manifold. We extract 
the large N behavior of the space-time extent in our model and compare the result with 
the prediction obtained by the low energy effective theory. Another dynamical issue to 
be addressed in this context is the space-time uncertainty relation, which was proposed as 
a principle for constructing nonperturbative string theory pi |. We extract the large N 
behavior of the space-time uncertainty of our model and confirm that the model indeed 
satisfies the proposed principle. 

Another attractive feature of the IIB matrix model as a nonperturbative definition of 
string theory is that its only parameter g is a simple scale parameter 0. One has to tune 
g suitably as one sends iV to infinity, so that the correlation functions have finite large N 
limits. According to Ref. || , Wilson loop operators can be interpreted as the string creation 
and annihilation operators, and it was found that g 2 N should be fixed in order to obtain the 
light-cone string field Hamiltonian in the large N limit. It is a non-trivial test of the model 
to verify if the correlation functions of Wilson loops really have a universal large N scaling. 
We address this issue in the present model and show that there is indeed a universal large 
TV scaling at fixed g 2 N '. 

We also address yet another important dynamical issue in this model, namely the ques- 
tion of equivalence to ordinary super Yang-Mills theory in the sense of Eguchi and Kawai 
PJ, which is exactly the way large N reduced models first appeared in history. The cru- 

2 This means in particular that the string coupling constant, which is related to the vacuum expectation 
value of the dilaton field, is not a tunable parameter. We come back to this point in Section 4.3. 
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cial observation is that large N gauge theory does not depend on the volume (under some 
assumptions), which inspired Eguchi and Kawai to propose the zero- volume limit of large 
N gauge theory as a model equivalent to the gauge theories in an infinite volume 0. One 
of the assumptions is that the (1<n) d symmetry of the model is not spontaneously broken, 
where D is the space-time dimension. However, in the purely bosonic case in D > 2, the 
symmetry is spontaneously broken at weak coupling |2"2]| , thus preventing one from taking 
a continuum limit. This led to modifications of the model |2"2" , [23] , |2~3] , |2"5] , [Ztj| ] so that the 



(Ztv) d symmetry is not spontaneously broken while keeping the equivalence valid. In the 
super symmetric case, the effective action which induces the spontaneous symmetry breaking 
of the (1*n) D symmetry is naively cancelled by the contributions of fermions. Indeed, in the 
scalar field case, it has been shown that the reduced model is equivalent to the field theory 



without such modifications [p7| . We observe in the present supersymmetric model that the 
Eguchi-Kawai equivalence indeed holds at least in a finite range of scale. What is rather 
remarkable is that actually this is true also for the bosonic case, which is contrary to what 
has been expected. 

In Section we describe the model we are going to investigate. In Section |3| we study 
the space-time structure of the model. In Section [| we present our results for correlation 
functions of Wilson loop and Polyakov line operators, and we discuss the Eguchi-Kawai 
equivalence as well as the universal scaling behavior. Section [| is devoted to a summary and 
discussion. In Appendix [A] we comment on the algorithm we used for the simulation. In 
Appendix [B] we present the corresponding results for the bosonic case for comparison. 

2 The model 

The model we investigate is a supersymmetric matrix model obtained by dimensional reduc- 
tion of 4D SU(iV) super Yang-Mills theory. The partition function is given by 



Z = J dA e~* b J dipdip e 
Sb = — "r^tr[A M , A/] 2 , 



Sf = -^tr(^(n a/3 [A M ,^]), (2.1) 



3 



where (/x = 1, . . . , 4) are traceless N x N Hermitian matrices, and ip a , ip a (a = 1, 2) are 
traceless N x N complex matrices. The measure is defined as 



dip dip 



dA 



n 

a=l 
4 

n 



i,j=l i=l i=l 

N N 

n{dRe(^) JJ dIm(^) lJ }n{d(^)45(^(A 



(2.2) 

v — , (2 ' 3) 

_i<j i=l i=l 

This model is invariant under 4D Lorentz transformations 0, where A^ transforms as a vector 
and ip a as a Weyl spinor. are 2x2 matrices acting on the spinor indices, and they can 
be given explicitly as 

Ti = iax, T 2 = ia 2 , T 3 = ia 3 , T 4 = 1. (2.4) 
The model is manifestly supersymmetric, and it also has a SU(iV) symmetry 



Ay. — > VA^ ; ip a ^Vip a V^ ; ip a ^Vip a V\ 



(2.5) 



where V G SU(iV). All these symmetries are inherited from the super Yang- Mills theory be- 
fore dimensional reduction. The model can be regarded as the four- dimensional counterpart 
of the IIB matrix model [[J. 

In contrast to unitary matrix models, where the integration domain for the partition 
function is compact, the first nontrivial question to be addressed in Hermitian matrix models 
in general, is whether the model is well-defined as it stands. The problem can be most clearly 
understood by decomposing the Hermitian matrices into eigenvalues and angular variables, 
where a potential danger of divergence exists in the integration over the eigenvalues, even 
at finite N. This issue has been addressed numerically for the supersymmetric case [H| at 
iV = 3 as well as the bosonic case [12] up to iV = 6. Exact results are available for N = 2 



There is also a perturbative argument which is valid when all the eigenvalues are well 
separated from each other |L3|, This reasoning agrees with the conclusions obtained for 
small iV [12, 14, |2£|. In particular, supersymmetric models in D = 4,6, 10 are expected to 
be well-defined for arbitrary N. Our simulations confirm that this is indeed the case for 
D = 4. 



3 When one defines the IIB matrix model nonperturbatively, a Wick rotation to Euclidean signature in 
needed. This is also the case for the present model. Hence by Lorentz invariance we actually mean rotational 
invariance. 
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Since the model is well-defined without any cutoff, the parameter g, which is the only 
parameter of the model, can be absorbed by rescaling the variables, 

A, = g l/2 X^ (2.6) 
^ Q = g 3/ ^ a - (2.7) 

Therefore, g is a scale parameter rather than a coupling constant, i.e. the g dependence 
of physical quantities is completely determined on dimensional grounds. The parameter g 
should be tuned appropriately as one sends N to infinity, so that each correlation function 
of Wilson loops has a finite large N limit. Whether such a limit really exists or not is one 
of the dynamical issues we address in this work. 

We now discuss the Eguchi-Kawai equivalence J7|, which is the equivalence between 
reduced models and the corresponding gauge theories in the large N limit. In its proof 
based on the Schwinger- Dyson equation, one has to assume quantities of the type 

(tr(e tt " A ")tr(e-' fc " A '0) + 0) (2.8) 

to vanish. Assuming in addition large N factorization, the vanishing of fl2.8|) is equivalent 
to (tr(e ifcAlj4 '')) = 0, which is guaranteed if the eigenvalues of are uniformly distributed 
on the whole real axis in the large N limit. The fact that the present model is well-defined 
without any cutoff implies that the eigenvalue distribution of A^ is not uniform, but it has 
a finite extent for finite N. Hence, the Eguchi-Kawai equivalence is quite nontrivial even in 
the supersymmetric case. Here the situation is more subtle than in the case of the unitary 
matrix model version of a large N reduced model J7j. There, the model has the (%n) D 
symmetry — > e 2mm ^/ N u^ ( m ^ = 0, 1, • ■ ■ , AT — 1), hence quantities like (tr(U^) n ) vanish, 
unless the symmetry is spontaneously broken. 

One might be tempted to consider a model defined by the partition function ( |2.1| ) but 
without imposing the traceless condition on A^. We denote such a model as the U(iV) model, 
to be distinguished from the original model, which we call the SU(iV) model. The U(iV) 
model has the U(l) 4 symmetry 

A M -> A^ + c^Itv, (2.9) 

where a M is a real vector. Note, however, that the trace part of A^ in the U(iV) model simply 
decouples because A^ appears in the action only through commutators. The transformation 
fl2.9|) acts on the decoupled trace part and hence it cannot play any physical role. Indeed 
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the quantity Q2.8p calculated with the U(iV) model or with the SU(iV) model is exactly the 
same. Thus, considering the U(7V) model does not help. 

Next we comment on the method we use to study the model. Details can be found in 
Appendix A. The integration over fermionic variables can be done explicitly and the result 
is given by det A4, A4 being a 2(N 2 — 1) x 2{N 2 — 1) complex matrix which depends on 
A^. Hence the system we want to simulate can be written in terms of bosonic variables as 

dA e~ Sb detM . (2.10) 

A crucial point for the present work is that the determinant det Ai is actually real positive, 
as we prove in Appendix A. Due to this property, we can introduce a 2(N 2 — 1) x 2(N 2 — 1) 



Hermitian positive matrix D = M'M, so that det M. = Vdet V, and the effective action of 
the system takes the form 

S cS = S b -UndetV . (2.11) 



We apply the Hybrid R algorithm |L8| to simulate this system. In the framework of this 
algorithm, each update of a configuration is made by solving a Hamiltonian equation for a 
fixed "time" r. The algorithm is plagued by a systematic error due to the discretization of r 
that we used to solve the equation numerically. We performed simulations at three different 
values of the time step At. Except in Fig. |2|, we find that the results do not depend much 
on At (below a certain threshold), so we just present the results for the value At = 0.002, 
which appears to be sufficiently small. 

We also note that there is an exact result 

(trF 2 ) = -(tr(^[A„^] 2 )) = 6g 2 (N 2 - 1), (2.12) 

which can be obtained by a scaling argument, similar to the one used for the bosonic case 
T3fl . We used this exact result to check the code and the numerical accuracy. 



3 The space-time structure 

We first study the space-time structure of the reduced model. In the IIB matrix model, 
the eigenvalues of the bosonic matrices A^ are interpreted as the space-time coordinates 



19, 20]. However, since the matrices A„ are not simultaneously diagonizable in general, 



the space-time is not classical. In order to extract the space-time structure, we first define 
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the space-time uncertainty A by 

A2 = ~ ^m^ {{UA " U%r ' (31) 

i 

which is invariant under Lorentz transformation and SU(iV) transformation (|2.5|) [O. This 



formula has been derived in Ref. |13|] based on analogy to quantum mechanics, regarding 
as an operator acting on a space of states. It has the natural property that A 2 vanishes 
if and only if the matrices are diagonalizable simultaneously. For each configuration 
A^ generated by a Monte Carlo simulation, we maximize EdiUA^Wju} 2 with respect to 
the SU(iV) matrix U . We denote the matrix which yields the maximum as t^ max , and we 
define x^ = {U^^A^U^^u as the space-time coordinates of iV points {i = 1, • • • , N) in 
four- dimensional space-time. 

Note that x^ should be identified with the dynamical variables denoted by the same x^ 
in Ref. [19]. There, the bosonic matrices A^ and the fermionic matrices ip a are decomposed 



into diagonal and off-diagonal elements as 

(*P a )ij = iaAj + Paij ( ( Paii = 0)- (3.2) 

The off-diagonal parts a^j and <p a ij are integrated out using the "Lorentz gauge" in the one- 
loop approximation, which is valid when the points x^ (i — 1, • - • , N) are well separated from 
each other. Thus one obtains the effective action for x^ and which can be considered as a 
low-energy effective action of the supersymmetric large N reduced model. In order to get the 
effective action only for x^, one still has to integrate over which cannot be done exactly 
for D = 6 and D = 10 (IIB matrix model). In D = 4, however, the integration over £ Qi can 
be carried out exactly and the system of x^ is described by a simple branched polymer with 
an attractive potential between the points connected by a bond. In D = 6 and D = 10, 
the system of x^ is expected to be described by some complicated branched-polymer like 
structure. Thus, although the one-loop approximation might seem quite drastic, the low 
energy effective theory of x^ still has a nontrivial dynamics. In Ref. |l||, some plausible 



mechanisms for the collapse of the x^ distribution in IIB matrix model have been discussed. 
What we have described in the previous paragraph provides a way to extract the low-energy 
effective theory of x^ from the full model without perturbative expansions. In particular, we 
can check explicitly whether the one-loop approximation adopted in Ref. ]T9| really captures 
the low energy dynamics of the supersymmetric large N reduced model. 
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We first look at the distribution p(r) of the distances r, where the distance between 
two arbitrary points x,- L 7^ Xj is measured by a/ (x~i — Xj) 2 . In Fig. |] we plot the results 



for iV = 16, 24, 32 and 48. We first note that the distribution at small r falls off rapidly 
below r/y/g ~ 1.5, independently of N. (This behavior is also seen in the bosonic case 
shown in Fig. [H].) This observation is in agreement with the argument in Ref. |19| that the 
ultraviolet behavior of the space-time structure of the model is controlled by the SU(2) matrix 
model. There, this argument has been used to justify the introduction of a ^-independent 
ultraviolet cutoff in the low energy effective theory, which otherwise suffers from ultraviolet 
divergence due to coinciding x^s. Our observation confirms that the ultraviolet cutoff 
is indeed generated dynamically if one treats the full model nonperturbatively instead of 
making perturbative expansions around diagonal matrices. 

P(r) 
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0.008 
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Figure 1: The distribution of distances p(r), plotted against r/^/g for N = 16, 24, 32 and 
48. 

In both, the supersymmetric as well as the bosonic case, we observe that the distribution 
shifts towards larger r as one increases N. In order to quantify this behavior, we define the 
extent of space-time by 




We denote this quantity by R new in order to distinguish it from the definition of the extent 
of the space-time R = 



(j^ti(A^)) used in Ref. fl3| . R, which roughly corresponds to 

J (J °° dr r 2 p(r)), is logarithmically divergent in the 4D supersymmetric case due to the 
asymptotic behavior p(r) ~ r~ 3 at large r [|15]. On the other hand, -R new does not suffer 
from this divergence as eq. ( |3.3|) shows. In Fig. ^| we plot the results for the space-time 
extent -R new as well as those for the space-time uncertainty y/ (A 2 ) for N = 16, 24, 32 and 
48. We repeat the same measurements for the bosonic model with N up to 256 and include 
the results in Fig. |2| for comparison. We see that the effect of fermions enhances R ne w and 
suppresses y/ (A 2 ) considerably. However, the power of the large iV behavior does not seem 
to be affected. 
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Figure 2: R ne w/y/9 and a/ (A 2 )/g, plotted against AT. The results for the bosonic model 
are also included. The lines are fits to the power behavior oc iV 1//4 , which is predicted 
theoretically. (In the labels we use a short-hand notation.) 



Let us discuss the results for -R ne w In the bosonic case, the data can be nicely fitted to a 
power behavior with R nC w/ \J~9 = 1.56(1) -iV 1 ' 4 . As expected, the observed large A" behavior of 
.Knew is the same as the one obtained for R in Ref. [|13|]. In the supersymmetric case, the large 
A" behavior of -R new can be predicted by the branched polymer picture based on the one-loop 
approximation Since the Hausdorff dimension of branched polymers is four, g?h = 4, 

the number of points A" grows as the extent -R new of the branched polymer, A" ~ (R new / £) dH . 
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Here i is the minimum length of the bond, which is of O(^fg) as we have already discussed. 
Thus one obtains i? n ew ~ y/gN 1 / 4 . The data in Fig. |2| seem to be consistent with this 
prediction. Fitting the data to this power behavior, we obtain R n ew/-\/g = 3.30(1) ■ iV 1 / 4 . 

One might be surprised that supersymmetry does not affect the power of the large N 
behavior of the space-time extent R aew . We recall, however, that in the bosonic case the 
explanation is completely different — although the power is the same ||13|| . There the one- loop 
perturbative expansion around diagonal matrices yields a logarithmic attractive potential 
between all the pairs of eigenvalues. The one-loop effective potential is dominant as far as 
the extent of the eigenvalue distribution is larger than ^fg N l / A . One can therefore put an 
upper bound on the space-time extent R < ^fgN 1 ^. What happens actually is that this 
upper bound is saturated. The behavior R ~ y/~g N 1 ^ can also be shown to all orders in the 



1/D expansion (13] 



Let us turn to the results for the space-time uncertainty. Using the one-loop perturbative 
expansion, (A 2 ) can be roughly estimated as 

The powers of -R new and a/ (A 2 ), as well as the coefficients we observe, are in qualitative 
agreement with this estimation. The bosonic case has been studied before in Ref. |13[. The 



data in Fig. can be nicely fitted to a power behavior with ^(A 2 )/^ = 0.907(3) • iV 1 / 4 . 
Thus, in the bosonic case we obtain a/ (A 2 ) ~ 0.58i? now |13| , which indicates a signficant 



deviation from the classical space-time picture. On the other hand, in the supersymmetric 
case we obtain ^/(A 2 )/g = 0.730(2) • iV 1//4 , hence our result amounts to a/(A 2 ) ~ 0.22 R ncw , 
coming closer to the classical space-time picture. 

We will see in the next section that the scale parameter g should be taken to be 0(1/ y/~N) 
in order to obtain a universal scaling behavior for the Wilson loop correlators. This means 
that the space-time uncertainty in the physical scale remains finite, rather than vanishing, in 
the large N limit. Therefore the present model satisfies the space-time uncertainty principle 



proposed for nonperturbative definitions of string theories [^T 

4 Wilson loop correlation functions 

In the interpretation of a large N reduced model as a string theory, Wilson loop operators 
correspond to string creation operators ||. Therefore, the existence of a non-trivial large 
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N limit of the Wilson loop correlators is an absolutely crucial issue. It has been addressed 
before in the 2D Eguchi-Kawai model, where non-trivial large N scaling has indeed been 
observed [l( 

We define the 'Wilson loop" and the "Polyakov line" operators as 

W(k) = l t r(e ifcX V fcX2 e- ifcXl e- ifcX2 ) , 

P(k) = ^tr(e^), (4.1) 



where are dimensionless matrices defined in eq. (2.6). For convenience we have chosen 
particular components of in the above definitions, but the choice of the directions becomes 
irrelevant when taking the vacuum expectation value, due to Lorentz symmetry and parity 
invariance. In the actual calculations we take an average over all possible choices of the 
components in order to enhance the statistics. 

The real parameter k represents the dimensionless "momentum" that characterizes the 
momentum density distributed along the string. The physical (dimensionful) momentum 
variable is given by k phys — kj ^/g. We have to tune g depending on N, so that the correlation 
functions of the above operators have definite large N limits as functions of k p ^ ys . In the 
following, we always set g — 1 for N = 48 without loss of generality. 

In all plots except for Fig. |], we further assume g to be proportional to 1/y/N. This 
turns out to be consistent with large N scaling, hence g oc 1/ \/N can be regarded as one of 
our observations. 

4.1 One-point function and Eguchi-Kawai equivalence 

In this subsection we discuss the one-point functions, and we start with the Wilson loop 



(W(k)}. Also Ref. |16| presents some recent results about this quantity. 
In the small k regime it can be expanded as 

(W(k)) = l + ^ 4 (tr([X 1 ,X 2 ] 2 )) + 0(A; 6 ) 



1 - ifc 4 f iV - i- ) +0(k 6 ), (4.2) 



4 V N 



where we have used the exact result ( |2.12| ). Therefore, in order to make the small k regime 
scale, we have to take g oc 1/y/N, as we mentioned above. In Fig. |^ we plot (W(k)) against 
k/^fg. The small k region scales as it should, and the results agree with the analytical 



prediction (4J2). Moreover the scaling extends up to k/^/g = 0(1). 
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Figure 3: The Wilson 1-point function (W), plotted against A; p h y s = k/^fg. 
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Figure 4: The Wilson 1-point function (W) is plotted now logarithmically against k 2 /g, in 
order to visualize the extent of the area law behavior. The scale parameter g has been tuned 
as described in the text. 
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If the model is equivalent to ordinary gauge theory — namely to 4D pure super Yang-Mills 
theory with four supercharges — which is confining, then the Wilson loop should exhibit an 
area law behavior. In order to illustrate this behavior, we show a logarithmic plot of (W(k)) 
versus the area k 2 /g in Fig. |4j. In this figure only, we fine-tune g as a function of iV so that 
the scaling in the intermediate regime of k becomes even better. We stay with the convention 
0(48) = 1 and use the optimal values g{32) = 1.291, g{24) = 1.563, 0(16) = 1.929, which 
is not far from g oc 1/ y/N. The small deviation can be understood as a manifestation of 
finite N effects. Fig. |] shows indeed a region of k that corresponds to the area law behavior 
(W(k)) ~ exp(— const. k 2 ). Surprisingly, the area law behavior is also observed in the bosonic 
model, as Fig. [12] shows, which is quite contrary to what one might have expected [13|, |16 



In both cases, supersymmetric and bosonic, it is not clear from the data whether the area 
law extends to k = oo in the large N limit. We will discuss the observed area law behavior 
from a theoretical point of view later. 
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Figure 5: The Polyakov 1-point function (P), plotted against k/yfg. 



We now proceed to the one-point function of the Polyakov line. In the 2D Eguchi-Kawai 
model [l(J this quantity vanishes due to (Zjv) d symmetry. In the present model, however, 
there is no exact symmetry that could make such a quantity vanish, as we explained in 
Section 0. Note for instance that P(k = 0) = 1 for any configuration. Fig. |5] shows the 
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results for (P(k)). It falls off rapidly as k increases Q. Again we observe a good scaling with 
g oc 1/vJV. Remember also that (P(k)) is actually just a Fourier transform of the eigenvalue 
distribution. Therefore, the value of k at which (P(k)) drops to zero, which we denote as 
k , should be inversely proportional to the space-time extent R new . The observed scaling 
with g oc 1/yN is consistent with our result in the previous section, -R new ~ y/gN 1 ^. The 
result for the bosonic case is shown in Fig. [T3|. We obtain a similar behavior except for 
some oscillations in the large k region. In particular, scaling is confirmed with g oc 1/y/N. 
The value of ko is larger than the supersymmetric ko, as expected. The ratio of ko in the 
two cases is indeed roughly the inverse of the corresponding ratio of R ncw (the bosonic k is 
about twice as large as the supersymmetric one). 

The above observations concerning (P(k)) and -R n ew have an interesting implication on 
the Eguchi-Kawai equivalence. We recall that from the results for (W(k)), we phenomeno- 
logically concluded that the Eguchi-Kawai equivalence holds at least in a finite range of scale 
for both, the supersymmetric and the bosonic case. We would like to understand this from 
a theoretical point of view. As we mentioned in Section |2], in the proof of Eguchi-Kawai 
equivalence, (P(k)) is assumed to vanish. We have found that (P(k)) is indeed small for 
k > k , but not for k < k . This means that the proof works for k > k Q , but not for 
small k, which corresponds to the ultraviolet regime in the corresponding gauge theory. We 
also observed that ko remains finite with respect to a physical scale in the large N limit. 
A complementary understanding can be obtained by taking Gross-Kitazawa's point of view 
J24|. As explained in Ref. []13| |, the extent of the eigenvalue distribution of determines 



the momentum cutoff of the corresponding gauge theory p4 |. The observation in Section [3] 
that -R n ew ~ \fg ' N 1 ^ implies that the momentum cutoff remains finite in physical scale as 
iV — > oo. Let us assume that the momentum cutoff is finite, but large enough to attract the 
renormalization flow to the fixed point which corresponds to the universality class of gauge 
theory. Then the flow follows closely the renormalization trajectory of the gauge theory, in 
a certain regime. That would explain why the equivalence holds at least in a finite range of 
scale. However, since the momentum cutoff does not go to infinity in the large N limit, it 
is conceivable that the renormalization flow will leave the renormalization trajectory of the 
gauge theory at some low-energy scale eventually. In this case the observed area law would 
not extend to k = oo even in the large N limit. 



4 One may consider the small k expansion here, as in eq. (4.2). The result is (P(k)} = 1— 2^/c 2 (tr(X 1 2 )) + 
• • •. The fact that (tr(A At 2 )) is logarithmically divergent means that actually (P(k)) has a non-analytic 
behavior ~ 1 + const, k 2 In Ifcl around k = 0. 
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4.2 Multi-point functions and universal scaling 

In this subsection we proceed to the large N scaling of multi-point functions of Wilson loops. 
We first note that in the bosonic case, there are analytical results to all order in the l/D 
expansion [[Uj. The statement is that 

(0102 ■ ■ • O n ) con ~ O {^j^—^j for the bosonic case, (4.3) 

where Oi denotes a Wilson loop or a Polyakov line as defined in eq. (|4.1| ), and (• • -) con means 
that only the connected part is taken. The correlation functions should be considered as 
functions of /c p h ys = k/^/g, where g is taken to be proportional to \/\fN. Our results for 



the bosonic model shown in Figs. [TT] to clearly confirm this analytical prediction. Let us 
consider a wave-function renormalization for each operator, of en ^ = ZOi, so that connected 
correlation functions of the renormalized operators of en ^ become finite in the large N limit. 
Relation ( |4.3|) means, however, that we cannot make all the multi-point functions finite. If 
we make the two-point functions finite by choosing Z ~ O(N), then all the higher-point 
functions vanish in the large N limit. In the supersymmetric case, we will see that scaling 
is observed again with g oc 1/y/N, but in contrast to the bosonic case a universal Z that 
makes all the correlators finite seems to exist. In the following, we always set Z(N = 48) = 1, 
without loss of generality. 

Let us start with the two-point functions, for which we measure the following two corre- 
lation functions, 

Gf\k) = ({ImW^)} 2 > 

G { 2 P) (k) = ({ImP(k)Y). (4.4) 

We take the imaginary part in order to avoid subtraction of a disconnected part. [] (Note in 
this regard that since ImW(k) and ImP(k) are parity odd, the one-point functions (ImW(k)) 
and (ImP(A;)) vanish due to parity invariance of the model.) The results are shown in Figs. 
H and |7], respectively. If we multiply the data by (iV/48) 2 , they scale nicely with g oc 1/y/N. 
As a three-point function, we measure 

Gf\k) = ({lmW{k)) 2 ReW{k)) - ({ImW {k)) 2 ) (ReW {k)) . (4.5) 



We also measured a number of multi-point functions, which are not presented here since the relative 
errors are rather large. 
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ure 6: The Wilson 2-point function G 2 , multiplied by Z 2 oc N 2 , plotted against kj \f§- 




We multiply the data either by (iV/48) 3 , which is required for the universal scaling of all 
the multi-point correlation functions, or by (iV/48) 4 , which is the factor predicted for the 
bosonic model. The results are compared in Fig. ||. We do observe a nice scaling behavior 
with a factor of (iV/48) 3 , but the scaling becomes worse for a factor of (iV/48) 4 . 




s 



e-06 



e-08 



e-12 




Figure 8: The Wilson 3-point function multiplied by Z 3 oc A^ 3 , plotted against kj ^fg 

on the left. On the right we show the corresponding plot using the bosonic prediction 
Z 3 oc iV 4 instead, which leads to an inferior level of scaling. 



Similarly, four-point function we measure 

Gf\k) = ((ImW(k)) 4 } -3((lmW(k)) 2 }. 



(4.6) 



We multiply the data either by (iV/48) 4 , which is required for the universal scaling of all 
the multi-point correlation functions, or by (iV/48) 6 , which is the factor predicted for the 
bosonic model. The results are compared in Fig. |9|. Again the scaling behavior obtained 
with the factor for universal scaling is superior over the behavior with the bosonic factor. 
To summarize our results concerning Wilson loop correlators, we observe that 



(0)~O(1) 
(O l 2 ---O n ) l 



for n > 2. 



(4.7) 



These correlators scale as functions of /c p h ys = k/^/g, where g is taken to be proportional 
to 1/ yN . This means that all the multi-point functions of the renormalized operators 



O 



(ren) 



ZOi become finite in the large limit if we set Z ~ O(N), in contrast to the 



bosonic case. We will discuss further the implications of this universal scaling behavior in 
the next subsection. 
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Figure 9: The Wilson 4-point function G^' 1 , multiplied by Z 4 oc iV 4 , plotted against k/y/g 
on the left. On the right we show the corresponding plot using the bosonic prediction 
Z A oc iV 6 instead, which leads to an inferior level of scaling. 

Finally we comment on large the N factorization. In ordinary gauge theory, large N 
factorization can be shown by weak-coupling expansion as well as strong-coupling expansion. 
In a large N reduced model with Hermitian matrices, one cannot do a weak-coupling or a 
strong-coupling expansion, because g is not a coupling constant but a scale parameter, as 
we have mentioned. Hence large TV factorization is nontrivial. In the bosonic case, large iV 
factorization holds to all orders of the 1/D expansion ||13|| . Our observation (|4.7|) implies 



(0 1 2 ...O n ) = {0 1 ){0 2 )...{O n ) + 0[ — 



where the O (1/iV 2 ) contributions are due to (Oi0 2 ) conies) ■ ■ ■ (O n ), etc. Therefore the large 
N factorization is also valid in the supersymmetric case. 

4.3 Interpretation of the large iV scaling 

In this subsection, we further discuss the physical significance of the large N scaling (|4.7|) 
we observed. 

If one views the supersymmetric matrix model considered here as a non-perturbative 
definition of type IIB string theory, and the Wilson loops as fundamental strings, string 
unitarity requires a large iV behavior of the form N aXn for the connected correlators of n 
Wilson loops, where Xn = 2 — n is the Euler characteristic of the worldsheet. In order to 
compare our results for the supersymmetric case ( f4.7|) as well as for the bosonic case ( |4.3| ) 
to this behavior, we first drop the extra 1/N n factor in ( |4.3| ) and (¥77), which is due to the 
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chosen normalization ( |4.1[ ) of the operators 0j. Then we find that the connected correlators 
of Wilson loops change from an 0(N Xn ) behavior to an 0(1) behavior by the introduction 
of supersymmetry. Our results for the supersymmetric case indicate a = 0, which is not in 
contradiction to the above requirement of string unitarity, but it is an extreme case where one 
is far away from a perturbative expansion in genus. It would be interesting to understand the 
result a = analytically, directly from the matrix model. From the string theoretical point 
of view this indicates that the supersymmetric matrix model might automatically realize a 
kind of double scaling limit, as it has been known from (ordinary) matrix models of two- 
dimensional quantum gravity. In these models contributions from all genera are important. 

While we do not presently have an analytic understanding of the difference between the 
connected correlators in the supersymmetric and the bosonic case, we can try to make an 
educated guess, based on the perturbative expansion fl3.2|). 

Let us consider the bosonic case. After integrating out the off-diagonal elements pertur- 
batively, schematically each diagram gives a contribution 




where Xi denotes the diagonal elements as usual, and F, L, V3, V4 are the number of index 
loops (faces), propagators (links), 3-point and 4-point vertices, respectively. They obey the 
relations 

F + V -L = x , 41/ 4 + 3^3 = 2L , V = V 3 + V A . (4.10) 

Here \ is the Euler characteristic of the diagram given as % = 2 — 2h — n, where h is the 
genus (the number of handles in the diagram) and n is the number of the operators. 

Let us now integrate over the diagonal elements Xi, under the assumption that the infrared 
region dominates. This assumption implies that the x^s are not close. In fact we assume 
that they are generically separated by some scale R, which is a measure of the extension of 
our universe. Viewing the Xi's as the coordinates of the worldsheet in target space, the above 
assumption implies that the worldsheet is rough: points are scattered quite randomly. We 
can now estimate the value for the diagram considered by simply replacing the integration 
over the x^s by their typical separation, assuming the existence of a suitable measure which 
implements the above hypothesis. We obtain 




(4.11) 
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where we have used R ~ yjg N 1 ^ . So the assumption of infrared dominance reproduces the 
observed large N behavior. 

In the supersymmetric case, the fermion diagonal elements make the power-counting more 
complicated. However, let us naively consider the contribution ( |4.9|) . When we integrate 
over the diagonal elements Xi, let us assume that the ultraviolet rather than the infrared 
region dominates. Namely, we assume that smooth worldsheets are favoured dynamically 
due to supersymmetry Q More precisely we assume that the dominant Xi configurations are 
such that Xj's, which are connected in a given diagram, are as close together as the actual 
distribution of Xj's allows. Above we have seen that there seems to be such a minimal length 
£, which is of the order ^/g, characterizing the distribution of \x{ — Xj\. Replacing \xi — Xj\ 
with this minimal length in (|4.9|) we arrive at 

Thus the assumption of ultraviolet dominance leads to our observed results in the super- 
symmetric case. Of course it remains to be understood why there is a difference between 
infrared and ultraviolet dominance in the bosonic and supersymmetric cases. 

If the above naive argument is true, it also implies that the contributions in the super- 
symmetric case do not depend on the genus h either, and consequently that diagrams of 
all topologies contribute with equal weight. This is reminiscent of the double scaling limit 
of matrix models. For example, let us consider a Hermitian one-matrix model with the 
partition function 

Z = Jdcj) e"^ 2 "^ 3 ), (4.13) 

being a N x N Hermitian matrix. The n-point correlation functions of loop operators 
behave as 

where e is the parameter related to the k and to the coupling constant A as 

k~\ , A-A c ~e 2 , (4.15) 

and A c is the critical coupling constant. Note that the large N behavior for fixed e is iV x , 
which we obtain for the bosonic large iV reduced model. When one takes the large iV limit 



6 It should be remarked here that such an increased smoothness of the bosonic part in a supersymmetric 
worldsheet has been observed in toy models for superstrings pm . 
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with Ne>' 2 — g st l fixed, the dependence on h disappears and the remaining power behavior 
e~ n can be absorbed into the wavefunction renormalization of each operator. In the large 
N reduced model, we do not have the parameter corresponding to A. The large iV scaling 
we observed for the supersymmetric case is formally the same as one obtains in the double 
scaling limit of the one-matrix model. In this sense, one might say that the double scaling 
limit is taken automatically in the supersymmetric large N reduced model, and that the 
string coupling constant g stT is not a tunable parameter but is fixed dynamically. If this is 
the case, it should be considered as a very satisfactory feature of the supersymmetric large 
N reduced model as a nonperturbative definition of superstring theory, since we do expect 
the string coupling constant g stT , which is related to the vacuum expectation value of the 
dilaton field, to be fixed dynamically (if superstring theory is treated nonperturbatively). It 
is also interesting that the qualitative difference of the large N behavior for the bosonic and 
the supersymmetric case might be traced back to the smoothness of the worldsheet, which 
is itself a dynamical question to be addressed. This point needs further clarification. 

5 Summary and discussion 

In this paper, we have studied the large N dynamics of a supersymmetric large N reduced 
model by means of Monte Carlo simulations. 

We studied the space-time structure represented by the eigenvalues of the bosonic ma- 
trices. In particular, we found that the large N power behavior of the space-time extent is 
consistent with the branched-polymer picture based on the one-loop perturbative expansion 
around diagonal matrices. The effect of fermions in the space-time extent was observed 
by the enhancement of the coefficient in the power behavior, but not in the power itself. 
The power appears to be the same for the bosonic and supersymmetric case. We empha- 
sized, however, that the theoretical explanation is completely different. We also found that 
the space-time uncertainty is clearly reduced for the supersymmtric case, which means that 
space-time comes closer to the classical behavior. Even in the supersymmetric case, the 
space-time uncertainty is found to be finite in the physical scale in the large N limit. We 
argued that this implies that the model satisfies the uncertainty principle for the nonpertur- 
bative definition of superstring theory. 

The large N scaling behavior of Wilson loop correlators is observed at fixed g 2 N. Al- 
though this scaling of g is the same as in the bosonic model, there is a striking difference from 
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the bosonic case in the wave-function renormalization with the multi-point functions. In the 
bosonic case, there was no universal scaling behavior: keeping two-point functions finite, all 
the higher-point functions vanish. In the supersymmetric case, we observed a clear trend for 
all the higher-point functions to become finite in the large N limit. We gave a perturbative 
argument that this result for the supersymmetric case might be understood if we assume 
smooth worldsheets to dominate. This argument also implies that all the topologies of the 
worldsheet contribute with equal weight to the amplitude. All these features are reminiscent 
of the double scaling limit of matrix models. 

We also addressed the issue of Eguchi-Kawai equivalence. By searching for the area law 
behavior in the one-point function of the Wilson loop, we concluded that the equivalence 
does hold at least in a finite region of scale. What is rather surprising is that the area 
law behavior has been observed also for the bosonic model. This suggests that the bosonic 
model is also equivalent to ordinary large N Yang-Mills theory at least in a finite region of 
scale, which is contrary to what has been generally believed. We argued, however, that this 
conclusion can be understood from a more theoretical point of view based on the large N 
behavior obtained for R 11CW and the one-point function of the Polyakov line. It is an open 
question whether this equivalence extends to the far infrared regime. 

To summarize, we have gained new insight into the dynamical properties of the large 
iV behavior of a supersymmetric large N reduced model. We hope that our findings shed 
light on the dynamical aspects of the most interesting 10D version of our model, i.e. the 
IIB matrix model. In this respect, it is encouraging that the large N scaling of Wilson 
loop correlators in the present model has been observed at fixed g 2 N, which coincides with 
the result obtained by requiring that the loop equations of the IIB matrix model should 
reproduce the string field Hamiltonian. We presume that a large iV scaling of Wilson loop 
correlators — like the one we observed — also holds for the IIB matrix model; then the 
only difference would be the spontaneous breakdown of Lorentz symmetry. One of the good 
news revealed in the present work is that low energy effective theory, based on the one-loop 
approximation, does already capture the low energy dynamics of the supersymmetric matrix 
model. We therefore hope to address the most interesting issue of spontaneous breakdown 
of Lorentz invariance by using the low energy effective theory — which is in 10D far more 
complicated than in the 4D case. We are going to report on Monte Carlo studies of IIB 
matrix model along these lines in the near future. 
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A The algorithm for the Monte Carlo simulation 

In this appendix, we explain the algorithm we use for the Monte Carlo simulation of the 
super symmetric matrix model. Only in this appendix we set g = 1 for simplicity. 

We first carry out the integration over fermionic matrices to obtain the explicit formula 
for the fermion determinant. We calculate 

Z f [A] = J di()d$ e~ s f, (A.l) 

where we use the notation introduced in eq. (|2.1|). We define a set of generators t a e gl(iV,C) 
by 

(*")« = Mm. (a = l,---,N 2 ), (A.2) 
where i a and j a are integers running from 1 to N, specified uniquely by 

a = N(i a - 1) + j a . (A.3) 

We also introduce the notation a = N(j a — l) + i a . The fermionic matrix ip a can be expanded 
in terms of t a as 

N 2 

faUij = J2*Paa (A.4) 
a=l 

where ip aa = {ip a )i a j a - ^Pa and can be expanded similarly with the coefficients ip aa = 
(4>a)i a j a and A aix = (A^) iaja . Note also that A^ = {A afl )* due to the Hermiticity of A^. 
We define the structure constants g a b c of gl(iV,C) by 

g abc = ti(f[t a ,t b }) 

= fijaib^jbicfijcia ~ ^jdb^jbia^jaic- (A-5) 
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The fermionic action then reads 



Sf — ~g a bc4 , ca(X^)a/3Aafj,1pb/3 

= -TpaaM'^plphp, (A.6) 

where 

M' aa w = -gaUT^^A^. (A.7) 

We first integrate out (ipa) nn and (4>o)nn using the 5 functions in the measure (|2.2|). We 
get a factor of 1/N followed by the replacements 

N-l N-l 

(4>a)NN => ~ ^2(?Pa)jj ] {^o)nN ~ ^2$a)jj (A. 8) 

3=1 3=1 

in the fermionic action. The integration over the remaining Grassmann variables yields 
det M., where M is the 2(A^ 2 — 1) x 2(A^ 2 — 1) complex matrix defined by 

M aa)b(3 = M' aa>b/3 - M' N 2 a ^5 iaja - M' aatN2 p5 ibjb (A.9) 

(the indices a and b run from 1 to A^ 2 — 1). Thus, we obtain 

Z f [A] = ± det M. (A. 10) 

We first want to show that the determinant det M. is real positive []. For this purpose we 
note that the matrix M. satisfies the identity 02M.02 — A4*. Hence if (p aa is an eigenvector of 
Jvt with an eigenvalue A, then ip aa = {^ap^ap)* is an eigenvector of M. with an eigenvalue 
A*. It is important that the two vectors ip aa , ip aa are linearly independent. The determinant, 
which is the product of all the eigenvalues of Ai, should therefore be real and positive semi- 
definite. In the case of 6D or 10D (IIB matrix model) versions of the supersymmetric 
large A" reduced model, the fermion integral yields a complex effective action in general. 
This causes the notorious sign problem, which makes standard Monte Carlo simulations 
practically inapplicable for large N. In the present case, since the determinant det M. is real 
positive, we can introduce a 2(A^ 2 — 1) x 2(A^ 2 — 1) Hermitian matrix D = Ai^Ai, which 
has real positive eigenvalues, and det Ai = V det T>. Therefore we have written the effective 
action for the bosonic matrices A^ in eq. Q2.ll ) as 



S cS = S b -UndetV . (A.ll) 



7 This has been already reported in Ref. p3 as a numerical observation. For related work, see Ref. |E9| 
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We apply the Hybrid R algorithm JTB| to simulate this system []. 

The first step of the Hybrid R algorithm is to apply the molecular dynamics method [[31 
We introduce a conjugate momentum for A a ^ as X ail , which satisfies X a ^ = (X afJj )*. The 
partition function can be re-written as 



dXdA e~ H , (A.12) 
where H is the "Hamiltonian" defined by 

H = \Y, x *n x *» + S bi A \ - \ IndetP. (A.13) 



2 

fia 



The update of X afJ , can be done by just generating X aAt with the probability distribution 
exp(— \ ^2 \X a ^\ 2 ). In order to update A atl , we use the Hamiltonian equations 



dA afl (r) dH 



dr dX afl 



X- a », (A.14) 



dX a Jr) dH 1 / dV ^ A <9S fe _ 

- -tr — — V~ x - — (A.15) 



(A.16) 



dr dA afl 2 J dA a ^ 

Along the "classical trajectory" given by the Hamiltonian equation, 

(i) H is invariant, 

(ii) the motion is reversible, 
(hi) the phase-volume is preserved, 

d(A(r),X(r)) _ 
d(A(0),X(0)) 

where (A(t),X(t)) is a point on the trajectory after evolution from (A(0),X(0)). Therefore, 

generating a new set of (A, X) by solving the Hamiltonian equation for a fixed "time" interval 

r satisfies detailed balance. This procedure — together with the proceeding generation of 

X atl with the Gaussian distribution — is called "one trajectory", which corresponds to "one 

sweep" in ordinary Monte Carlo simulations. 

In order to solve the Hamiltonian equation numerically, we have to discretize the "time" r. 

A discretization which maintains the properties (ii) and (iii) is known. The slight violation of 

(i) for finite At causes systematic errors. One can in principle eliminate the systematic error 

completely, by making a Metropolis accept/reject decision at the end of each trajectory. But 

in the present case, the overhead for this procedure is rather large. We therefore decided to 

8 Ref. [[[o) gives an overview of effective algorithms for dynamical fermions, including the Hybrid R algo- 
rithm. 
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omit that step, and just use a sufficiently small At. Still we can use the specific discretization 



of Ref. [ 18| ) which we explain below, to minimize the systematic error. As we explain later, 



we do find a good convergence in small At, and the systematic error is well under control. 
We introduce a short-hand notation for the discretized X a/J (r) and A a ^{r), 

X§=X a „{rAr) ; A<$ = A afl (sAr) . (A.17) 

The Hamiltonian equations are discretized as 

A) _ i(0] ,^v(0) 
liafM — /i a/Jj -f- s^an 

(n+ i) (n+ i) At {n) 

4 (™+i) , 3At {m) 

^a/i — ^ 2 

= x £ } + A -{^ +l) -^(^ +l) )} ' (A - 18) 

where n = 0, 1, • • • , i/ — 1, m = 1, • • • , v — 1, and -R^ t +2 ' ) is defined by 

^ = $ aa I ^bf3 , (A.19) 

d(4^ } )$ = A^t (A (;+i) )r? . (A .2o) 



Here r\ aa are complex variables generated by the Gaussian distribution exp(— J2 aa \Vaa\ 2 ) 
The judicious choice of the argument of JVO is the tool to reduce the systematic error [0 
We solve eq. ( A. 20 ) with respect to $ by means of the conjugate gradient method [32 



which is iterative. Each iteration involves a multiplication of the matrix T> with some vector 
v. Since T> is a 2(N 2 — 1) x 2(A^ 2 — 1) matrix, storing T> requires 0(N 4 ) memory, and 
multiplying T> with v naively involves 0(N 4 ) arithmetic operations. Actually we can do 
much better than this. We first recall that V = M^M, where M. is the 2(N 2 — 1) x 
2(N 2 — 1) matrix defined in eq. (|A.9|) . The point is that the number of nonzero elements of 



M. is only 0(N 3 ) (not 0(A^ 4 )). Indeed, the multiplication Aiv can be done economically 
as follows. 

We consider 

W aa = M aab pV b/3 , (A.21) 
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and define the quantities w' aa and v' aa , where a runs from 1 to N 2 as in Ai', by 



v' aa = v aa for a = 1, ■ ■ ■ , N 2 — 1 , 



v'n*« = -][>««, ( A ' 22 ) 

1a = 3a 

<a = M' aab pv' b p. (A.23) 



Now w nn can be written as 



w aa ={ W ™- W N*« for u ia r Ja (A.24) 
1 w aa otherwise. 



Thus the problem reduces to calculating the matrix- vector product in eq. (|A.23|) . Using 
definition ( |A.7|) , we obtain 

Wij = F' t )ap[A»v , () ] j i, (A.25) 

where w' a and v' a are N x N matrices associated with w' aa and v' aa , respectively, as in eq. 
( |A.4| ). The commutator in eq. ( |A.25| ) requires 0(N 3 ) arithmetic operations. Thus we save 
O(N) operations. In addition, we do not have to store neither g abc nor A4. Multiplication 
of MO with some vector v is done in the same way. 

A similar technique should be used to calculate R afl in eq. ( |A.19|) . Note first that it can 
be written as R Cfl = T Cfl + (T Sfl )*, where T CAt is given by 



- * act \ rj A J ) 

C A» / aab/3 



\dA c 

^aa = (M aab ^ bf3 y . (A.26) 



We define $' and 1 J r ' in terms of $ and as we defined v' in terms of v before in eq. (|A.22 ). 



Now we can re-write T CAt as 



T cfl = gj- (< a (M') aabf} & b p) . (A.27) 



Using again eq. ( |A.7| ), we obtain 

(TX = -(r,) Q/3 [<,$^ , (A.28) 

where &' a , ^>' a and are N x N matrices associated with & aa , ^' aa and T' respectively, 
as in eq. (|A.4j). 



There are two parameters v and At in this algorithm. We can choose z/Ar so that a 
typical autocorrelation time is minimized. We have taken z/Ar = 1 throughout the present 
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work, and v = 200, 280, 500 for each of the cases N = 16, 24, 32, and v = 500, 600 for 
iV = 48. Except in Fig. |2|, we observed that the results are reasonably well converged 
at v — 500, At = 0.002, so we just present those results. For Fig. [| we carried out an 
extrapolation to At = by assuming the At dependence of some observables Q(At) to be 

Q(At) - Q(At = 0) ~ (At) 2 ■ (tr(A 2 )) Ar . (A.29) 

This assumption has been checked for (trF 2 ) with the exact result ( [2.12| ). We also observed 
that (tr(A 2 ))& T behaves as 

(tr{A 2 )) Ar ~ ci - c 2 log At , (A.30) 

for small At, where C\ and c 2 are constants depending on N. [] This implies that it diverges 
logarithmically for At — > 0, which is consistent with the theoretical prediction discussed 
below eq. ( |3.3| ). 

Let us comment on the required computational effort of our algorithm. The dominant 



part comes from solving the linear system (|A.20|) using the conjugate gradient method. First 



of all, we find that the number of iterations necessary for the convergence of the method 
seems to grow linearly with the size of the matrix V, namely as 0(N 2 ). This is much worse 
than the full QCD case with a fixed quark mass, where the number of iterations does not 
depend on the system size. We may interpret this phenomenon as a sort of "critical slowing 
down", since the present system corresponds to QCD in the chiral limit. As we have seen, 
the number of arithmetic operations for each iteration is of order N 3 . Therefore, the required 
computational effort of our algorithm is estimated to be 0(N 5 ). 



For the bosonic case, we use the heat bath algorithm in the way proposed in Ref. [[13 
which requires an effort of 0(iV 4 ). We note, however, that application of a Hybrid Monte 
Carlo algorithm |33|] allows for an 0(N 3 ) algorithm for the bosonic case, which might be 
useful for proceeding to much larger N. 

Finally, we comment on the numbers of configurations used for the measurements. For 
the super symmetric case, they are 3060, 1508, 1296, 436 for iV = 16,24,32,48, respectively. 
For the bosonic case, we used 1000 configurations for each N. 



9 In QCD the At dependence of the systematic error is 0(Ar 2 ) A similar argument leads to the 



be (At) 2 log At. 



assumption (A.29). Due to eq. (A.30), the At dependence of the systematic error in our case is expected to 
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B Results for the bosonic case 



For comparison we show in this appendix the results for the bosonic case. By the bosonic 
case we mean a model obtained by just dropping the fermions from the supersymmetric 
matrix model described by the partition function ( |2.1| ). Fig. |iy shows the distribution p{r) 
defined in Section [| Figs. [11| to [T7| show the Wilson loop and Polyakov line correlators 
defined in Section ffl. We take g oc 1/ y/N (g — 1 for N = 48) and plot the results against 
^phys = k/y/g, as in the supersymmetric case. We multiply the results by (N / '48) 2(n_1 ) for 
n-point functions. 

The data scale nicely in agreement with the theoretical prediction for large N given by eq. 
( [4.3| ). For comparison we also show the 3-point and the 4-point Wilson loop correlators with 
the renormalization factors, which were used successfully in Sec. 4 for the supersymmetric 
case. We see very clearly that the bosonic prediction is the correct one in this case. 
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Figure 10: The bosonic distribution of distances p(r), plotted against rj^fg for N = 16, 24, 
32 and 48. 
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Figure 11: The bosonic Wilson 1-point function (W), plotted against k/y/g. In this case, 
the small k prediction amounts to 1 — (N/6)k 4 . 
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Figure 12: The bosonic Wilson 1-point function (W) is plotted now logarithmically against 
k 2 /g, in order to visualize the extent of the area law behavior. 
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Figure 13: The bosonic Polyakov 1-point function (P), plotted against k/y/g. 
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